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The octopus arm is a muscular hydrostat and due to its deformable and highly flexible 
structure it is capable of a rich repertoire of motor behaviors. Its motor control system uses 
planning principles and control strategies unique to muscular hydrostats. We previously 
reconstructed a data set of octopus arm movements from records of natural movements 
using a sequence of 3D curves describing the virtual backbone of arm configurations. Here 
we describe a novel representation of octopus arm movements in which a movement is 
charactehzed by a pair of surfaces that represent the curvature and torsion values of points 
along the arm as a function of time. This representation allowed us to explore whether the 
movements are built up of elementary kinematic units by decomposing each surface into a 
weighted combination of 2D Gaussian functions. The resulting Gaussian functions can be 
considered as motion primitives at the kinematic level of octopus arm movements. These 
can be used to examine underlying principles of movement generation. Here we used 
combination of such kinematic primitives to decompose different octopus arm movements 
and characterize several movement prototypes according to their composition. The 
representation and methodology can be applied to the movement of any organ which 
can be modeled by means of a continuous 3D curve. 
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INTRODUCTION 

Octopuses are considered to be among the most developed and 
intelligent animals in the invertebrate kingdom, where at least 
part of their skills can be attributed to the high maneuverability 
of their arms and the capacity of the peripheral nervous system 
to process sensory information and control arm movements. The 
octopus uses its arms for various tasks such as locomotion, food 
gathering, hunting, and sophisticated object manipulation (Wells 
and Wells, 1957; Fiorito et al, 1990; Mather, 1998). The versa- 
tile and adaptive nature of octopus movements is mainly due 
to the flexible nature of the octopus arms which do not con- 
tain any rigid skeleton. The octopus arm is a muscular hydrostat 
built of closely packed arrays of muscle fibers organized in three 
main muscle groups: parallel, perpendicular, and helical that runs 
obliquely to the long axis (Matzner et al., 2000). A constant vol- 
ume constraint that holds for muscular hydrostats allows forces to 
be transferred between the longitudinal and the transverse mus- 
cle groups. The movements of a muscular hydrostat are based on 
combinations of four elementary movements that can occur at 
any location: elongation, shortening, torsion, and bending (Kier 
and Smith, 1985). Therefore, both structural support and force 
transmission are achieved through the arm's musculature, such 
that the biomechanical principles governing octopus arm move- 
ments differ from those operating in arms with a rigid skeletal 
support. 



The octopus nervous system is divided into a central and 
peripheral nervous systems. Axial nerve cords are projecting from 
the brain along the center of each arm, and the peripheral neurons 
located in the axial nerve cords are organized into an exten- 
sive nervous system comprising both sensory and motor circuits 
(Young, 1971). Behavioral studies suggest that the nerve cord 
circuitry and the peripheral components play a major role in 
the control of the complex actions performed by octopus arms 
(Altman, 1971; Wells, 1978). 

Analyses of octopus reaching (Gutfreund et al., 1996, 1998; 
Sumbre et al., 2001; Yekutieli et al., 2005a,b) and fetching 
movements (Sumbre et al., 2001, 2005, 2006) have revealed 
some control principles that underlie movement generation. 
During reaching a bend point propagates along the arm fol- 
lowing an invariant velocity profile. Fetching movements use a 
vertebrate-like strategy, reconfiguring the arm into a stiffened 
qMfl5!-articulated structure. These movements were studied by 
analyzing the kinematics of the movements of specific points 
along the arm which display several stereotypical characteristics. 
Electromyographic recordings and detailed biomechanical simu- 
lations assisted in revealing common principles which reduce the 
complexity associated with the control of these movements. The 
travelling bend, used in arm extension movements, was found to 
be associated with a propagating wave of muscular activations, 
where simple adjustments of the excitation levels at the initial 
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Stages of the movement can set the velocity profile of the whole 
movement. Recently, a soft robotic arm inspired by the octopus 
arm has been designed in order to reproduce the octopus tentacle 
motor performance and to examine the possibility for the imple- 
mentation of motor control principles identified in the octopus 
as part of its control (Laschi et al, 2009; Calisti et al., 2011). 

However, describing the movements of specific points along 
the arm is insufficient for analyzing the fuU complexity of octo- 
pus arm movements. To determine whether the kinematics of 
octopus arm movements can be described by a reduced set of 
motion primitives requires analysis of different types of arm 
movements and of the shape of the entire arm as it moves through 
space. Motion primitives can be regarded as a minimal set of 
movements, which can be combined in many different ways 
giving rise to the richness of vertebrate and invertebrate move- 
ment repertoires and allowing motor learning of new skills (Flash 
and Hochner, 2005; Bizzi et al., 2008). Motor primitives have 
been inferred at various levels of the motor control system. Sub- 
movements were shown to be combined at the kinematic level 
(Krebs et al., 1999; Rohrer et al., 2002), a reduced set of static 
force field underlie controlling arm posture (Mussa-Ivaldi and 
Bizzi, 2000; d'Avella et al, 2003, 2006), while movement dynamics 
can be learned through a flexible combination of dynamic prim- 
itives (Thoroughman and Shadmehr, 2000). Dynamical move- 
ment primitives were also used to model attractor behaviors of 
autonomous non-linear dynamical systems and rhythmic move- 
ments (Ijspeert et al., 2002, 2013), and discrete and rhythmic 
movement elements were used to investigate single-joint and 
multi-joint motor behaviors (Sternad et al., 2000; Sternad and 
Dean, 2003). Inferring motion primitives from octopus arm 
movements may help understand underlying principles and kine- 
matic optimal measures, and provide new understanding of how 
the nervous system in muscle hydrostats handles the complexi- 
ties associated with the control of hyper-redundant arms. This 
may also facilitate designing control systems for hyper-redundant 
robotic manipulators. 

Here we refer to the behavioral level and aim at describing 
octopus arm behaviors as being composed of elementary kine- 
matic units to which we also refer to as motion primitives (Flash 
and Hochner, 2005). We believe that identifying basic kinematic 
patterns is the first step in further investigating the existence of 
primitives also at the control, movement dynamics, and muscle 
activation levels as well as the neural control levels as was demon- 
strated in earlier studies of the octopus motor system (Gutfreund 
et al., 1996, 1998; Sumbre et al, 2001, 2005, 2006; YekutieK et al., 
2005a,b). As was discussed in Flash and Hochner (2005) elemen- 
tary building blocks may exist at all the above levels of motor 
representation but the most immediate and direct way is to search 
for elementary units at the kinematic level. Movement strokes 
with specific spatial and temporal features and submovements 
were shown to successfully describe both periodic and discrete 
motions and were indicated as plausible building blocks of human 
and monkey movements (Sosnik et al, 2004; Polyakov et al., 
2009). Furthermore, in robotics research locomotion trajecto- 
ries for a humanoid robot were constructed based on kinematic 
motion primitives derived from humans' locomotion trajecto- 
ries (Moro et al., 2011, 2012). Relations between behavioral 



and control levels were suggested in different earlier studies, for 
example: hand trajectories of stroke patients were shown to be 
composed of submovements with velocity primitives obeying the 
minimum jerk model (Flash and Hogan, 1985) whose number 
was found to decrease as the patients gained better control of their 
limb (Rohrer et al., 2004), simple curved two-dimensional trajec- 
tories that follow the two-third power law (Lacquaniti et al., 1983) 
were described by means of parabolic units and corresponded to 
neural activation states identified using a hidden Markov Model 
(Polyakov et al., 2009). Another example consists of grasping and 
object manipulation movements described as arising from well- 
coordinated combinations of basic motor actions-arm transfer 
and hand shaping (leannerod, 1994). 

An algorithm for 3D tracking and analysis of octopus arm 
movements (Yekutieli et al, 2007; Zelman et al., 2009) enabled 
us to create a large database of many types of modeled octopus 
arm movements. Here we describe a new framework allowing 
the extraction of kinematic units from these reconstructions. 
The octopus arm movements were represented by a pair of sur- 
faces describing the curvature and torsion values of the arm. 2D 
Gaussians for each surface were extracted such that each Gaussian 
represented the characteristic shape of the curvature or torsion 
along a section of the octopus arm during some time interval. We 
found that Gaussian functions generally fit quite well the contin- 
uous form of the configurations that the octopus arm can take 
with respect to both the time and the arm index dimensions: the 
curvature and torsion values were observed to change smoothly 
along the arm length for any quasi-static arm configuration, and 
the magnitude of curvature or of torsion at any specific point was 
gradually changing with time during the movement. Gaussian- 
like functions were previously used in composing hand velocity 
(Thoroughman and Shadmehr, 2000) and limb position profile 
(Hwang etal, 2003). 

The resulting Gaussians were divided into clusters whose cen- 
troids defined kinematic units, and each movement was repre- 
sented as a weighted combination of such units. These kinematic 
units can be used to form a language of motion primitives, allow- 
ing characterization and representation of a large repertoire of 
octopus arm movements. We show how these kinematic units 
can be used to classify octopus arm movements into meaningful 
groups. Understanding how kinematic primitives can be utilized 
and combined can greatly contribute not only to studies of motor 
control in octopus arms and other hyper-redundant appendages 
but can also provide a deeper understanding of motor control 
systems in general. 

METHODS 

The analyzed octopus arm movements in our study were per- 
formed by four specimens of Octopus vulgaris, weighing 200 
(female), 200 (male), 450 (male), and 470 (female) g. The animals 
were maintained in 50 x 40 x 40 cm tanks containing artificial 
seawater. The water was circulated continuously in a closed sys- 
tem through a biological filter of Orion, gravel and coral dust. 
Water temperature was held at 16°C in a 12 h light/dark cycle. 
Prior to the video recording sessions, the animals were moved to 
a bigger glass tank (80 x 80 x 60 cm) with a water temperature 
of 18°C. 
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SPATIO-TEMPORAL REPRESENTATION OF MOVEMENT AS A PAIR OF 
CURVATURE AND TORSION SURFACES 

Since the octopus arm displays no well-defined landmarks, a 
skeletal representation can be naturally used to model the octopus 
arm using curves which prescribe its virtual backbone. The back- 
bone was found using a "grass fire" algorithm that extracts the 
middle line of the arm: first the contour of the arm is separated 
into two sides, dorsal and ventral, from base to tip. Then two dis- 
tinct waves are initiated from the two sides of the contour and are 
propagated at an equal speed inward. The set of points where the 
wave fronts collide is the midline (Yekutieli et al., 2007). 

The reconstructions of octopus arm movements result in a 
sequence of 3D curves prescribing the virtual backbone of the 
octopus arm as its configuration changes during the movement. 
Figure 1 presents an extension movement as a sequence of 60 3D 
curves that prescribe the virtual backbone during the move- 
ment. For each curve, green, and red points mark the base of 
the arm (that was aligned between sequential images using tex- 
tural cues) and the tip, respectively. Given a sequence of m 
3D curves as an input, we wished to construct a pair of sur- 
faces describing the values of the curvature and torsion along 
these curves. Since arm configurations were reconstructed from 
video records whose sample rate was 50 frames/s, the smoothness 
of the motion between consecutive configurations of a move- 
ment was guaranteed, and a spline function was used to smooth 
noisy points as necessary (Yekutieli et al., 2007; Zelman et al., 
2009). 

Each 3D curve was first represented by (n = 100) sample 
points uniformly distributed along the curve. This 3D curve was 
then approximated by a cubic smoothing spline constructed of 



piecewise third-order polynomials passing through the n sample 
points. Approximation was achieved, by considering both the 
smoothness of the spline and the distance between the spline and 
the sample points. Formally, given the data site x(j) and the cor- 
responding data values y(j) for j = !,...,«, the cubic smoothing 
spline / minimizes: 

where the integral over the second derivative of/ is over the small- 
est interval containing all the entries of x. The smoothing param- 
eter p defines the tradeoff between the success in approximating 
the data points and the smoothness terms. 

We then calculated the curvature (k) and torsion (t) values for 
the n sample points along each of the 3D curves. The curvature 
was calculated using the circle passing through three successive 
points as an approximation of the osculating circle to the curve 
at the middle point. This is formally described by Calabi et al. 
(1998): Let A,B,C be three successive points on the curve C such 
that the Euclidean distances are a = d(A,B), b = d{B,C), c = 
d{A,C). Let A denote the area of the triangle whose vertices are 
A,B,C, and let s = j(a + b + c) denote its semi-perimeter, so that 
A = -y/s(s — a)(s — b)(s — c). Then the radius of the circle pass- 
ing through the points A,B,C is computed leading to the formula 
for its curvature: 

v/s(s — fl)(s — b)(s — c) 

k(A, B, C) = 4— -. 

abc 

In this study we will use the word curvature as the inverse of the 
radius of curvature. 

The torsion along a 3D curve, defined as t = ^, was cal- 
culated for a pair of successive points as the angle between the 
normals to the planes defined by the successive triangles corre- 
sponding to these points, divided by the distance between the 
points (Boutin, 2000). Let A,B,C,D,E be five successive points on 
the curve such that habc and hcoE are the normals to the planes 
defined hj A,B,C and C,D,E respectively, and the Euclidean dis- 
tance between points B and D is d[B,D). Then the torsion at point 
C is calculated as: 

COS^^ {flABC-ncDE) 

^^^^ = d^M • 

Finally, the curvature and torsion values calculated for a sequence 
of 3D curves were represented on two surfaces that separately 
described the curvature and torsion as a function of time and arm 
index. The result was a pair of smooth and normalized curvature 
and torsion surfaces, such that a single arm movement was com- 
pactly represented by a pair of « by n matrices. This representation 
was invariant to rotation and translation in a Cartesian coordinate 
system, as curvature and torsion measures are intrinsic (i.e., they 
do not depend on the orientation and position of the arm in 3D 
space). 

Figure! presents the curvature and torsion surfaces for the 
extension movement presented in Figure 1. The relatively high 




X 



FIGURE 1 I A spatio-temporal profile of an extension movement 
shown as a sequence of the virtual backbone of quasi-static arm 
configurations. Tine presented 3D curves are the result of the 
reconstruction process in which the virtual backbone prescribing the 
octopus arm is detected (see section spatio-temporal representation of 
movement as a pair of curvature and torsion surfaces). The virtual backbone 
was found by a "grass-fire" algorithm, green and red points mark the base 
of the arm (that was aligned between sequential images using textural 
cues) and the tip respectively (Yekutieli et al., 2007). 
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Curvature surtace tor arm 



Torsion surtace tor arm 





FIGURE 2 I Curvature and torsion surfaces extracted for the movement shown in Figure 1. Tine values are given as a function of the arm index and time. 



values of the curvature surface generally describe propagation of 
a bending section from the middle of the arm toward the tip. The 
decrease in the torsion values as the movement proceeded means 
that the arm configuration became relatively planar. 

SURFACE DECOMPOSITION USING GMM 

Gaussian Mixture Model (GMM) is a statistical method for 
density estimation and data clustering (McLachlan and Peel, 

2000) . In this model a Gaussian fitting method can be used 
to approximate a function of one variable by a weighted 
sum of 1 -dimensional Gaussians. As GMM is a generalized 
framework, it can approximate any multidimensional data by 
a set of multivariate Gaussians. The model uses an itera- 
tive process which optimizes the Gaussians' parameters by 
the Expectation Maximization (EM) algorithm (Xuan et al., 

2001) . 

A function of one variable (y = f(x)) can usually be approx- 
imated by a mixture of ID Gaussians, where each Gaussian is 
defined by its mean and standard deviation. In our case, we refer 
to a surface as a function of two variables z =f(s,t), where z 
stands for either the curvature or torsion values, and s,t refer to 
the arm index and time dimensions, respectively. We therefore use 
2D GMM to approximate a surface by a weighted combination of 
2D Gaussians. Specifically, a surface is approximated as: 



zis, t) = ■ g[v.i,i:i](s, t). 



where ^ is a Gaussian defined by a 2 x 1 mean vector |i and 2x2 
covariance matrix E, and w is the Gaussian weight. The Gaussian 
g is defined as: 



2n 



172 ^^P 



(x-yiy^-\x-\i) 



where x ■■ 



The mean vector |x corresponds to the position 



the covariance matrix E correspond to the standard deviation 
of the 2D Gaussian. Its two eigenvectors correspond to the axes 
of the Gaussian with respect to a fixed coordinate system, such 
that the covariance matrix defines the shape and orientation of 
the Gaussian. 

We also added a criterion to choose the right number of 
Gaussians into which each surface should be optimally decom- 
posed, based on the Minimum Description Length (MDL) prin- 
ciple. The MDL descriptor to be minimized here is the Bayesian 
Information Criterion (BIG), as developed by Andrews and Lu 
(2001): 



BIG : 



-2 ■ L + d- log(n) 



of the Gaussian center on the surface, and the two eigenvalues of 



where L is the log likelihood of the mixture of Gaussians, d is 
the number of parameters in the model (number of degrees of 
freedom) and n is the number of observations in the sample. 
The BIG criterion allows choosing the most parsimonious model, 
i.e., the model which best describes the data with respect to the 
number of Gaussians it uses for the decomposition. [See also 
Bhat and Kumar (2010) for a more detailed derivation of the BIG 
formula] . 

Figures shows the approximation of a curvature surface 
by a weighted combination of four Gaussians. Intuitively, each 
Gaussian can be illustrated as a hill, whose center, shape, orien- 
tation, and height are defined by the Gaussian parameters. The 
decomposition into 2D Gaussians allows us not only to explore 
Gaussians as possible units enabling to define the kinematics 
of octopus arm movements, but also to compactly represent a 
surface as a weighted sum of 2D Gaussians. 

CLUSTERING ALGORITHM 

The GMM allows us to decompose octopus arm movements 
into curvature and torsion 2D Gaussians which describe their 
kinematics. We refer to each of the resulting Gaussians as a 
data point, whose dimension corresponds to the number of 
parameters defining a Gaussian (section surface decomposition 
using GMM). To cluster these points (i.e., the 2D Gaussians) 
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FIGURE 3 I Gaussian Mixture IVIodel of the curvature surface of an extension movement. (A) Tlie input surface. (B) The resulting mixture of Gaussians. 



into meaningful groups we used the kmeans clustering algo- 
rithm which is an unsupervised clustering method. The output 
of the kmeans algorithm is k disjoint clusters, where each cluster 
includes a different number of points and is represented by a cen- 
troid that can be regarded as the average of all the points assigned 
to the cluster, kmeans uses a two-phase iterative algorithm to 
yield a clustering result which minimizes the point-to-centroid 
distances summed over all k clusters. 

The distance usually employed for kmeans is the Euclidean 
distance (Hastie et al., 2009), but here we want to improve 
the clustering from two points of view. First, we designed a 
Weighted Euclidean Distance, i.e., for each of the parameters of 
the Gaussians — center, shape, area, and orientation, we separately 
computed the Euclidean distance among the different elements 
of the sample. We then got four distances, each being related to 
one of the four parameters. The quantity to be minimized in the 
kmeans algorithm at each step is then an average of these four dis- 
tances. Second, we used the Gap-Statistics (Tibshirani et al., 2001) 
as a criterion of the optimal number of clusters to be used. Gap- 
Statistics compares the within-clusters distance of the distribution 
(given by kmeans) to the within-distance W^^^ of a Monte-Carlo 
sample drawn within the range of the reference distribution. This 
criterion was used for example by Ben-Hur et al. (2002) and 
Pedersen and Kulkarni (2006). The idea of this approach is thus to 
compare the graph of log( Wjt) (log of the within -cluster distance) 
with its expectation under an appropriate null distribution. The 
mathematical rationale of this approach is explained in greater 
detail by Tibshirani et al. (2001). Defining B as the number of 
generated data sets, the Gap-Statistics is expressed as: 



Gap (k) = l/B log {^tb) - log(^*) 

b= 1 

The optimal number of clusters is the minimal k which gives a 
local maximum of the Gap. 



RESULTS 

Our data set consisted of 60 reconstructions of octopus arm 
movements that included extension movements. Extension in the 
octopus arm is generally characterized by a bend propagating 
along the arm (Figure 11). Some of these extension movements 
were preceded by initialization movements, referred to as pre- 
extension movements, in which the octopus arm moved from an 
initial random position to a configuration that seemed to be ideal 
for the extension (Figure 12). 

Careful examination of the video sequences allowed us to 
define start- and end-points of the extension phase as charac- 
terized in earlier studies (Gutfreund et al., 1996, 1998). In an 
extension movement, a bend is created usually near the base of 
the arm and is propagated along the arm toward the tip where 
it vanishes, while the base of the arm points in the direction of 
propagation. A pre-extension movement is generally defined as 
a movement in which an arbitrary configuration of the arm is 
reconfigured to an initial extension configuration. Based on these 
observations we initially divided our data into 25 pre-extension 
and 60 extension movements. In order to characterize sets of 
kinematic units and determine synthetic rules allowing recon- 
struction of the observed movements, we next decomposed the 
movements into curvature and torsion Gaussian units (as defined 
in section surface decomposition using GMM) and analyzed these 
units as described below. Since each movement was defined by 
a specific combination of kinematic units, we could classify the 
movements into sub-groups, such that all the movements in a 
sub-group were defined by a combination of similar kinematic 
units. In order to explain the different phases of the movement 
analysis as clearly as possible we focus here mainly on the group 
of extension movements. 

DECOMPOSITION AND CLUSTERING OF KINEMATIC UNITS 

Curvature and torsion surfaces were extracted for all the octo- 
pus arm movements in our database (see section spatio-temporal 
representation of movement as a pair of curvature and torsion 
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surfaces). The curvature and torsion values at the tip of the 
arm could be very high (and sometimes noisy) relative to the 
values along the proximal and middle parts of the arm and 
were, therefore, analyzed separately. The resulting surfaces were 
approximated using the GMM, yielding decompositions of each 
curvature/torsion surface into 2D Gaussians units (see section 
surface decomposition using GMM). These units were found to 
naturally describe the surfaces, as each unit essentially described 
bending or torsion along a defined section of the arm and its 
movement along the arm as function of time. 

A set of 2D Gaussians was assembled as a set of kinematic units 
separately for the pre-extension and extension movement groups. 
A set of Gaussians can be variously clustered by referring only 
to a subset of the parameters defining the 2D Gaussians, namely 
the center location, size, shape and orientation (see section clus- 
tering algorithm). These parameters were easily extracted from 
the mean and covariance matrix of a 2D Gaussian; the coordi- 
nates of the Gaussian center on the surface (i.e., the time point 
and the arm index at which the Gaussian reached its maximum) 
were directly defined by the Gaussian mean. The orientation of 
the Gaussian (the angle between the Gaussian axes and the axes 
of the fixed coordinate system) was defined by the eigenvectors of 
the Gaussian covariance matrix. The projection of the Gaussian 
on the plane was an ellipse whose size and eccentricity were 
defined by the eigenvalues of the covariance matrix and the ratio 
between them. Finally, the relative influence of the Gaussian in the 
decomposition in which it participated was defined by its weight. 
The clustering results presented here were obtained by referring 
to the Gaussian's center location (Gaussian mean), Gaussian's 



shape (ratio between the eigenvalues of the Gaussian's covariance 
matrix), and Gaussian's weight. 

Figure 4 presents the clustering results obtained for the cur- 
vature and torsion Gaussians of the 60 extension movements. 
Gaussians marked by the same color belong to a single cluster. 
Executing kmeans with the Gap-Statistics method (section clus- 
tering algorithm) resulted in three clusters both for the curvature 
and torsion Gaussians. Coordinates of the centroids of the various 
clusters are presented in Table 1. 

These results essentially suggest that all the Gaussians com- 
posing the curvature and torsion surfaces of the extension move- 
ments can be classified into three types according to the values 
of the Gaussian's center location and shape. Explicitly, the blue 
curvature cluster (Figure 4 left) represents curvature Gaussians 
defining curvature along the proximal section of the arm during 
the movement. Examining the orientation of these Gaussians as 
defined by the eigenvectors of their covariance matrices shows that 
there was a relatively small angle between each of the Gaussian 
axes and the axes of the arm-index — time coordinate system 
(Table 2), that is, the internal Gaussian axes almost aligned with 
the direction of the arm-index and time axes. This characteristic 
means that the section of the arm to which these Gaussians relate 
did not change during the movement; we therefore term them 
"fixed" Gaussians. We suggest that these Gaussians correspond to 
movements used to aim the base of the arm toward a target point 
during the extension movement (see discussion below). 

The green and magenta clusters represent curvature Gaussians 
that travel toward the tip of the arm during the extension and 
are probably associated with the main characteristics of the bend 




FIGURE 4 I Clustering results for the curvature (left) and torsion (right) Gaussians that were extracted from a group of 60 extension movements by 
GMM. Tine arm index and time coordinates of eacli cluster centroid are given in Table 1. 
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Table 1 | The arm index and time coordinates of the cluster centroids shown in Figure 4. 



Cluster no. 

1 (green) 

2 (nnagenta) 

3 (blue) 



Centroids of curvature Gaussians 

Arm index coordinate Time coordinate 



0.4627 
0.6647 
0.1355 



0.1574 
0.4502 
0.4543 



Cluster no. 

1 (green) 

2 (nnagenta) 

3 (blue) 



Centroids of torsion Gaussians 

Arm index coordinate Time coordinate 



0.2687 
0.6245 
0.4061 



0.2463 
0.3288 
0.7207 



Table 2 | The mean median (11.1/2) and standard deviation (a) for the orientation values (in degree) of the resulting curvature and torsion 
clusters shown in Figure 4. 



Orientation of curvature Gaussians 



Cluster no. 

1 (green) 

2 (magenta) 

3 (blue) 



45.8 
40.4 
19.9 



1^1/2 

28.8 

31 

11.2 



42.6 
35.9 
24.2 



Cluster no. 

1 (green) 

2 (magenta) 

3 (blue) 



Orientation of torsion Gaussians 



21.4 
32.2 
39.4 



1^1/2 

8.3 

22.3 

29.2 



31.1 
30.5 
35.7 



propagation in extension movements. The mean orientation of 
each of these clusters is significantly larger compare with the blue 
cluster that refers to the base of the arm (Table 2). A similar 
interpretation is valid for the resulting clusters of the torsion 
Gaussians, where the green and magenta clusters refer to torsion 
Gaussians in the early stage of a movement, and the blue clus- 
ter refers to torsion Gaussians at the end of movement. Table 2 
presents the median, mean and standard deviation of the orienta- 
tion values of the resulting clusters. These findings are supported 
by earlier analyses and simulations of the stereotypical character- 
istics of an extension movement (Gutfreund et al., 1998; Yekutieli 
et al, 2005b). 

SYNTHESIZING ARM BEHAVIORS FROM KINEMATIC UNITS 

The 2D Gaussians were clustered by the kmeans algorithm based 
on their mean vector and covariance matrix values. The centroid 
point of each cluster represented the center of the cluster, i.e., 
the point giving the minimum sum of distances from it to all the 
data points in that cluster. Since the mean vector and covariance 
matrix which define a Gaussian as a data point also apply to the 
centroid point, a centroid point essentially defines a representa- 
tive Gaussian for its cluster. Each of these Gaussians has a unique 
position, size and orientation, thus uniquely defining the octopus 
arm movement in 3D space. 

We therefore consider the curvature and torsion Gaussians 
defined by the resulting centroid points as kinematic units that 
can be used to generate a set of behaviors of the octopus arm. 
These are local time-dependent behaviors, since they refer to a 
specific section of the arm at a specific time during the move- 
ment. Figure 5 presents the three curvature Gaussians (left) and 
the three torsion Gaussians (right) defined by the centroid points 
of the clusters found for the extension movements (Figure 4). A 
curvature unit alone defines a planar arm behavior, as it defines a 
change in the curvature level along a section of the arm as a func- 
tion of time, with a zero value for the torsion associated with the 
arm. Coupling a curvature and a torsion unit, such that both of 
them refer to a common section of the arm, defines a 3D behavior. 



However, a torsion unit on its own is meaningless since apply- 
ing torsion on a straight line representing the backbone of the 
arm has no effect on its configuration. A torsion unit has no 
significant effect also when it is coupled with a curvature unit 
that refers to a different section of the arm. In general, «c curva- 
ture units and rij torsion units can define «c • "r 3D behaviors, 
and since the nc curvature units define «c planar behaviors where 
they are not coupled with any torsion unit, they can overall define 
«c ■ («r + 1) behaviors. Figure 6 presents some of the behaviors 
that can be defined by the curvature and torsion kinematic units 
extracted for the extension group (Figure 5). They are shown as 
sequences of quasi-static configurations in 3D space, where the 
red, black and blue curves represent the initial, intermediate and 
final configurations in the sequence, respectively. 

CLASSIFYING OCTOPUS ARM MOVEMENTS 

The kinematic units (curvature/torsion Gaussians) extracted for 
the extension and pre-extension movement groups can be used 
further to classify movements in a given group into different 
sub-groups according to the mixture of Gaussians composing 
their curvature and torsion surfaces. Intuitively, movements that 
were decomposed into weighted combinations of similar kine- 
matic units were classified into the same sub-group, as they were 
assumed to be characterized by a similar 3D behavior. 

We represented each of the movements in our data set by a 
weighted combination of the curvature and torsion kinematic 
units defined for the group of movements to which the movement 
belonged. That is, a movement m which was approximated by a 
pair of curvature (C) and torsion (T) surfaces C= X!> ' sf 
and T = wj ■ gj , was represented by a row vector of the 

weights: w= [{wp}, j^J^}]- We applied the kmeans algorithm 
(section clustering algorithm) to the set of vectors of weights 
corresponding to a group of movements, such that the input to 
the algorithm is a matrbc of weights, where the rows correspond 
to the analyzed movements and the columns to the Gaussians 
that were previously identified as curvature and torsion units. 
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time time 

FIGURE 5 I The curvature (left) and torsion (right) centroid Gaussians of the extension group of movements. Each Gaussian is essentiaiiy the centroid of 

one of the ciusters in Figure 4. 




FIGURE 6 1 Simulations of octopus arm behaviors in 3D space behaviors show the characteristics of extension movements — directing 

defined by the curvature and torsion kinematic units which the base toward a target, initiaiization, and propagation of the 

were extracted for the extension movement group. These bend. 



Each row practically defines a movement as a weighted sum of 
the elementary Gaussian units (Figure 7). The kmeans algorithm 
separated the movements into clusters, such that movements 
belonging to the same cluster shared a similar pattern of weights. 
That is, a cluster consists of movements that can be spanned by 
a similar weighted sum of the available curvature and torsion 
units. We therefore refer to each of the resulting clusters as a sub- 
group of movements that share similar characteristics of their 3D 
behavior. The centroid point of each of the resulting clusters was 
considered a representative pattern of a weighted combination of 



kinematic units, which defined the behavior of the sub-group of 
movements in 3D space. We refer to these different behaviors as 
movement prototypes. 

Each of the three pairs of curvature and torsion surfaces 
presented in Figure 8 is defined by a weighted combination of 
kinematic units, corresponding to one of the extension sub- 
groups we found by the process just described. A pair of curvature 
and torsion surfaces defined a prototype which characterized 
one of the sub-groups by simulating a sequence of 3D curves 
whose curvature and torsion values corresponded to the values 
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FIGURE 7 I The kinematic units can span the 60 extension movements. 

Eacln movement is expressed as a weighed sum of curvature (left) and 
torsion (right) units. Clustering tliese weiglits results with patterns of 
weights, such that each one defines a movement prototype in this group. 



given by the curvature and torsion surfaces (Figure 9). Although 
the curvature surface of the first and third prototypes (Figure 8 
top, bottom) share a similar topographic structure, the differ- 
ence between their curvature surfaces defines prototypes with 
different characteristics. Relative to the first prototype, the third 
prototype (Figure 9 right) describes an extension movement in 
which a higher level of torsion is observed along with the prop- 
agating bend, causing the 3D configuration to deviate from the 
movement plane. Furthermore, the higher weight of the proximal 
curvature Gaussian results in a higher level of curvature along the 
base of the arm. The values presented by the torsion surface for 
the first prototype decrease during the second half of the move- 
ment, meaning that the configuration of the arm tends to become 
more planar as the movement progresses. For the first proto- 
type (Figure 9 left) the arm section around the propagating bend, 
which takes higher curvature values, creates a loop during the ini- 
tial phase of the extension movement. Compared to the first and 
third prototypes, in the second prototype (Figure 9 middle) the 
propagating bend starts with a higher curvature value and lower 
curvature values for the proximal section of the arm. 

PRE-EXTENSION RESULTS 

The analysis described above was also applied to the pre-extension 
movement group. The movements in this group refer to the 
actions that the octopus arm was observed to perform just before 
the extension phase has started. The well-defined time point in 
which the bend starts to propagate along the arm has been used 
to define the time at which a pre-extension movement ends. The 
kinematic units extracted for this group and the arm behaviors 
showed some similarities but also some unique characteristics. 
Figure 10 presents the single prototype that was found for pre- 
extension movements. It appears to represent the initializing 
phase of the arm, in which the base is directed toward a target 



and the bend (which is propagated during the extension) is gen- 
erated. The initialization of the bend is achieved by generating 
movements corresponding to the curvature and torsion kine- 
matic units on the same mid-arm section. Such dynamics may 
be associated with a minimal loss of energy due to interactions 
with drag forces. Computer simulations of the movements using 
the dynamic model of the octopus arm (Yekutieli et al., 2005a,b) 
will help to farther explore and characterize this prototype with 
respect to muscle activation and energy expenditure. 

DISCUSSION 

By carefully watching the octopus arm movements in video 
sequences and identifying the time points bounding the extension 
phase, we were able to divide our data set of reconstructed arm 
movements into two main groups, pre-extension and extension 
movements. The analysis described here was applied separately to 
each of these groups but we have presented results mainly related 
to the extension group. Equivalent results for the pre-extension 
group are also available. 

Instead of the common representation of octopus arm move- 
ments in 3D Euclidean space, we modeled each arm movement 
using pairs of curvature and torsion surfaces. These surfaces 
essentially describe the curvature and torsion values at the sam- 
pled points along the virtual backbone of the octopus arm as 
a function of time and arm index. Such pairs of curvature and 
torsion surfaces provide a compact description of arm configura- 
tion which is independent of the arm location in 3D space and 
is invariant to rotation and translation. Most importantly, this 
approach can be used to demonstrate the existence of kinematic 
units or motor primitives in octopus arm movements. 

The characteristics of the surfaces led us to examine whether 
they can be meaningfully decomposed. We applied the GMM, 
suggesting the use of 2D Gaussians as building blocks approx- 
imating the curvature and torsion surfaces of a movement. 
These 2D Gaussians provided a mathematically quantified rep- 
resentation whose hilly shapes fitted well to the topographic 
characteristic of the surface. We thus have demonstrated a mean- 
ingful representation for octopus arm movements and a method 
GMM for decomposing the movements into well-defined build- 
ing blocks (2D Gaussians) allowing us to further examine them 
as possible kinematic units. We have also applied an alternative 
method which decomposes a surface to its fundamental surfaces 
by analyzing the principal curvature values at each point on the 
surface. These parameters allow defining eight fundamental sur- 
faces (e.g., peak, pit, ridge) that correspond to the topography 
of a surface (Yilmaz and Shah, 2005). Interestingly, we found 
the results extracted by this method to be very similar to those 
achieved using the GMM — the positions of the peak fundamen- 
tal surfaces were highly correlated with the mean values of the 
positions derived using the 2D Gaussians. 

Intuitively, being able to represent the arm by a 2D curvature 
Gaussian corresponds to the propagation of a bend point along 
a defined section of the arm and during a defined time inter- 
val. All the kinematic properties — the affected section of the arm, 
the time interval, the maximal curvature value and the velocity 
of propagation — are simply defined by the center location of the 
Gaussian, by its covariance matrix and by the weight assigned 
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FIGURE 8 I Each pair of curvature and torsion surfaces defines one of the three prototypes of movements into which the extension movements were 
classified. Eacli surface is essentially a weighted combination of the curvature/torsion Gaussians extracted earlier. 



to this Gaussian. By clustering Gaussians with similar character- 
istics, we were able to characterize each cluster by its centroid 
and use the Gaussian representing the entire cluster as a stereo- 
typical kinematic unit. We obtained a set of such curvature and 
torsion units for each of the pre-extension and extension move- 
ment groups. Curvature and torsion units were then combined 
to simulate new movements in 3D space and to examine whether 
the entire observed repertoire of complex 3D octopus arm behav- 
iors can be spanned using the derived basic set of kinematic units. 
We found that patterns of weighted combinations of the kine- 
matic units can be clustered into prototypes of movements in the 



pre-extension and extension phases, allowing classification of the 
movements into sub-groups. 

The combinations of kinematic units which define the proto- 
types needs further investigation to reveal the principles underly- 
ing the execution of the different arm movements. Sumbre et al. 
(2001) have suggested that a relation between kinematic features 
and basic motor programs (embedded within the neural circuitry 
of the arm itself) greatly simplify the motor control of the octo- 
pus arm. In addition, a simple command producing a wave of 
muscle activation in a dynamic model was sufficient to repli- 
cate the kinematic characteristics of natural reaching movements 
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Ethograrti: pattern 2 



Ethogram: pattern 3 




FIGURE 9 I Three prototypes represent the sub-groups into which the 60 extension movements were classified. These prototypes, defined by tliree 
pairs of curvature and torsion surfaces (Figure 8) sliow tlie differences in the various extension movements. See text for further explanation. 
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FIGURE 10 I One prototype was found to represent the pre-extension 
movements. Pre-extension movements are intuitively understood as an 
initialization phase, during which the bend propagating during the extension 
phase is generated and the base of the arm is directed toward the target. 
The form of this prototype suggests that the movement is initialized by 
generating a new bend at the appropriate position by propelling the 
mid-section of the arm. 



(Yekutieli et al., 2005a). Specifically, it was found that natural 
extension movements can be generated by a dynamic model, in 
which a simple propagating neural activation signal is sent to con- 
tract muscles along the arm. In the model, the control of only two 
parameters fully specified the extension movement: the amplitude 
of the activation signal and the activation travelling time, such 
that different levels of activations can result in desired kinematics 
(Yekutieli et al., 2005b). We suggest that values of these two con- 
trol parameters can be associated with the characteristics of the 
kinematic units extracted here. That is, the weight, shape, orien- 
tation, and size of a Gaussians can be related to the amplitude of 
the activation signal and the activation travelling time. 

The relation between the kinematic units to the biomechanics 
of the octopus arm has to be examined (Feinstein et al., 2011). 
The arm morphology points to the dorsal group of the longitudi- 
nal muscles being much thicker than the ventral group, and both 
groups differ from the lateral groups. This anisotropy suggests 



that while bending movements to the left and right directions 
might be similar, this is not the case when comparing between 
upward, downward and sideward directions. The oblique mus- 
cles are composed of three pairs of helical bands, such that the 
handedness of the helix of one member of the pair is opposite to 
that of the other member of the same pair (Kier and Stella, 2007). 
This isotropy with respect to the arm axis supports that torsion 
toward the two different directions is applied in a similar manner. 

The results from our analysis agree with our data on octopus 
arm movements. The extension movement shown in Figure llA 
matched with prototype no. 2 of the extension group (Figure 9 
middle), as a movement in which a highly curved bend along 
a relative short section of the arm propagated rapidly toward a 
target. The lower movement Figure IIB matched with prototype 
no. 1 of the extension group (Figure 9 left), as a movement in 
which the arm moved relatively slowly toward the target with 
relatively low curvature values for the propagating bend. The 
Gaussians referring to the main characteristic of extension move- 
ments strengthen the previous findings of Gutfreund et al. (1996) 
of a stereotypic profile of the position and velocity of the bend 
point. They suggested that the position of the bend in space and 
time is a controlled variable, which simplifies motor control. The 
travelling bend, associated with a propagating wave of muscle 
activation (Gutfreund et al., 1998), was simulated as a biome- 
chanical mechanism in a dynamic model of the octopus arm 
(YekutieU et al, 2005a). 

Two pre-extension movements shown in Figure 12 matched 
with the prototype of the pre-extension group (Figure 10). A 
substantial manipulation of the initial arm configuration was 
involved by creating a bend in the arm and directing it toward 
the target. In this movement, curvature and torsion kinematic 
units are both applied on the mid-section of the arm during the 
pre-extension phase. These results demonstrate that the Gaussian 
description of movement primitives allows us to describe a com- 
plex motor behavior. Clearly, additional types of octopus arm 
movements other than the pre-extension and extension move- 
ments analyzed here are part of the motor repertoire of octo- 
pus behavior. While reconstruction and analysis of these other 
movements will probably reveal additional kinematic building 
blocks, we can expect that the general characteristics of octopus 
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FIGURE 11 I Two extension movements. The upper movement (A) 
was classified as prototype no. 2 of tine extension group (Figure 9 
middle), as a movement in whicin a highly curved bend was rapidly 
propagated toward a target, while the base of the arm stayed 
oriented with a fixed direction. The lower movement (B) was 



classified as prototype no. 1 of the extension group (Figure 9 left), 
as a movement in which the bend showed lower curvature values 
and moved relatively slowly toward the target while the direction of 
the base of the arm was not preserved. The movements progress 
from left to right in each panel. 




FIGURE 12 I A pre-extension movement as a sequence progressing from 
the upper left (A) to the lower right (F) frames. A substantial manipulation, 
creating a bend and directing it toward a target, was applied to the initial 



configuration (upper left). This movement is matched with the prototype of 
the pre-extension group (Figure 10). Frame (F) presents a temporal 
configuration that matches the beginning of an extension movement. 
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FIGURE 13 I A curvature surface which also refers to the elongation during an extension movement (left). Tine values of the ratio between the lengths of 
the proximal section of the arm (from base to bend point) and the length of the entire arm as a function of time (right). 



arm movements, mainly the smoothness in which curvature and 
torsion values vary along both the arm and the time dimen- 
sions wiU hold for other types of movements. Therefore, we 
beHeve that Gaussian functions could be efficiendy used also in 
the decomposition and description of those movements. 

Our results fit with Yekutieli et al.'s (2007) observations on 
the kinematic characteristics of the initiation of a reaching move- 
ment. The kinematic description is sufficiently rich for describing 
complex arm movements, although factors such as the biome- 
chanics of the octopus arm (e.g., the different type of muscles 
and the constant volume constraint), water drag forces and energy 
expenditure also strongly influence the arm movement charac- 
teristics. For example, the perpendicular drag coefficient for an 
octopus arm is nearly 50 times larger than the tangential drag 
force coefficient. This most likely affects the preferred arm config- 
uration during extension movements; only a small part of the arm 
is oriented perpendicularly to the direction of movement, mini- 
mizing drag (Yekutieli et al., 2005a). Yekutieli et al. (2005b) also 
found that the control of extension movements can be specified 
by the amplitude of the muscle activation signals and the activa- 
tion travelling time. The primitives we suggest here can be used to 
further investigate the relation between the kinematic and muscle 
activation levels. 

In our analysis the curvature and torsion surfaces were 
extracted for arm configurations whose length has been normal- 
ized. Replotting the curvature and torsion surfaces while showing 
the actual arm length values as analyzed from live data (Figure 13 
left) shows that the proximal section of the arm elongates dur- 
ing an extension (i.e., the section between the base and the bend 
point). This is demonstrated clearly in Figure 13 (right) which 
shows the ratio between the length of the proximal section to the 
length of the entire arm during an extension movement. Arm 
elongation has recently been shown to play a key role in the 
biomechanics and control of octopus arm movements (Hanassy, 
2008). Modeling the travelling bend along extension movements 
based on the propagation of muscle activation and stiffening 
wave (Gutfreund et al., 1998), where co-contraction of both the 



longitudinal and transverse muscles pushes the bends forward 
(Yekutieli et al., 2005b). Therefore, different ratios between the 
activation levels of longitudinal to transverse muscle can be used 
to control the elongation of the arm along the proximal section 
between the base and bending point. Gaussian units, which were 
found in this study to describe the travelling bend during exten- 
sion movements, will be further examined in order to support 
recent findings related to the biomechanics and control of the 
octopus arm. 

Our analysis presents a possible language of kinematic 
primitives — 2D Gaussians of either curvature or torsion which 
define and classify octopus arm behaviors by their different com- 
binations. Constructing a taxonomy of possible movements for 
a species is one approach to the study of its behavior. To con- 
struct a taxonomy of octopus arm movements and to reveal 
how combinations of components result in a variety of behav- 
iors Mather (1998) used components which consist of movements 
of the arm itself, the ventral suckers and their stalks, as well 
as the relative position of the arms and the skin web between 
them. Comparing similar movement taxonomies and ethograms 
(catalog of body patterns and associated behaviors) in the squid 
and various octopus species (Hanlon et al., 1999; Huffard, 2007; 
Mather et al., 2010) suggests that behaviors may be conserved 
throughout the evolution of these species. Our results identify a 
number of kinematic units, possible time-dependent units and 
sub-groups (Table 3). As more reconstructed octopus arm move- 
ments become available (Yekutieli et al., 2007; Zelman et al., 
2009), we will better be able to use our analytical tools to define 
a comprehensive language of motor primitives that incorporates 
the underlying kinematic principles, thus enriching the ethogram 
and taxonomy of octopus arm behavior. 

Our analysis provides a new framework for research on the 
kinematics and control of any natural or mechanical flexible 
manipulator. Possible arm behaviors can be simulated by syn- 
thesizing new combinations of the extracted Gaussians. New 
primitives can be hypothesized and tested on dynamic mod- 
els of the octopus arm and the resulting movements can 
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Table 3 | The analysis applied to each movement group yielded a number of kinematic units defining possible local-temporal 3D behaviors of 
the arm. 



Group 


Number of kinematic units 


Number of possible behaviors 


Number of sub-group prototypes 




Curvature 


Torsion 






Pre-extension 


1 


3 


3 


1 


Extension 


3 


3 


6 


3 



Each group was divided into a number of sub-groups which were represented by different movement prototypes. 



be compared with live movements. This, in turn, may allow 
future studies of activation commands at the neural control level, 
which may then enable operation of a real flexible manipulator 
to perform specified goal-oriented tasks (Laschi et al., 2009, 2012; 
Calisti et al.,2011). 
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